#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cstdint>
#include <ipps.h>
#include "private/throw.hpp"
#include "rtseis/postProcessing/singleChannel/detrend.hpp"

using namespace RTSeis::PostProcessing::SingleChannel;

class DetrendParameters::DetrendParmsImpl
{
public:
    void clear()
    {
        precision_ = RTSeis::Precision::DOUBLE;
        mode_ = RTSeis::ProcessingMode::POST_PROCESSING;
        linit_ = true;
    } 
    RTSeis::Precision precision_ = RTSeis::Precision::DOUBLE;
    RTSeis::ProcessingMode mode_ = RTSeis::ProcessingMode::POST_PROCESSING;
    bool linit_ = true; // This module is always ready to roll
};

class Detrend::DetrendImpl
{
public:
    DetrendImpl() = default;
    /// Copy constructor
    DetrendImpl(const DetrendImpl &detrend)
    {
        *this = detrend;
    }
    /// Deep copy operator
    DetrendImpl& operator=(const DetrendImpl &detrend)
    {
        if (&detrend == this){return *this;}
        parms_ = detrend.parms_;
        b0_ = detrend.b0_;
        b1_ = detrend.b1_;
        linit_ = detrend.linit_;
        return *this;
    }
    /// Destructor
    ~DetrendImpl() = default;
    /// Resets the module
    void clear()
    {
        parms_.clear();
        b0_ = 0;
        b1_ = 0;
        linit_ = true; // This module is always ready to roll
    }
    /// Sets the parameters for the class
    int setParameters(const DetrendParameters &parameters)
    {
        clear();
        parms_ = parameters;
        return 0;
    }
    /// Removes the trend from data
    int apply(const int nx, const double x[], double y[])
    {
        b0_ = 0;
        b1_ = 0;
        if (nx < 2){return 0;}
        computeLinearRegressionCoeffs(nx, x);
        #pragma omp simd
        for (int i=0; i<nx; i++)
        {
            y[i] = x[i] - (b0_ + b1_*static_cast<double> (i));
        }
        return 0;
    } 
    /// Removes the trend from data
    int apply(const int nx, const float x[], float y[])
    {
        b0_ = 0;
        b1_ = 0;
        if (nx < 2){return 0;} 
        computeLinearRegressionCoeffs(nx, x); 
        auto b0 = static_cast<float> (b0_);
        auto b1 = static_cast<float> (b1_);
        #pragma omp simd
        for (int i=0; i<nx; i++)
        {
            y[i] = x[i] - (b0 + b1*static_cast<float> (i));
        }
        return 0;
    }
    /*!
     * @brief Computes the coefficients for a linear regression
     *          \f$ \hat{y}_i = b_0 + b_1 x_i \f$
     *        using IPP functions where
     *          \f$ b_1
     *           = \frac{\mathrm{Cov}(x,y)}{\mathrm{Var}(x,y)} \f$
     *        and 
     *          \f$ b_0 = \bar{y} - b_0 \bar{x} \f$.
     *        The code is modified from:
     *        https://software.intel.com/en-us/forums/intel-integrated-performance-primitives/topic/299457 
     *        where I have assumed the data is generated by linearly
     *        spaced samples.
     * @param[in] length  Length of array y.
     * @param[in] pSrcY   Points at evaluated at indices.
     */
    int computeLinearRegressionCoeffs(const int length, const double pSrcY[])
    {
        double cov_xy;
        double mean_x;
        double mean_y;
        double var_x;
        // Mean of x - analytic formula for evenly spaced samples starting
        // at indx 0. This is computed by simplifying Gauss's formula.
        auto len64 = static_cast<uint64_t> (length);
        mean_x = 0.5*static_cast<double> (len64 - 1);
        // Note, the numerator is the sum of consecutive squared numbers.
        // In addition we simplify.
        var_x = static_cast<double> ( ((len64 - 1))*(2*(len64 - 1) + 1) )/6.
              - mean_x*mean_x;
        ippsMean_64f(pSrcY, length, &mean_y);
        cov_xy = 0; 
        #pragma omp simd reduction(+:cov_xy)
        for (int i=0; i<length; i++)
        {
            cov_xy = cov_xy + static_cast<double> (i)*pSrcY[i];
        }
        // This is computed by expanding (x_i - bar(x))*(y_i - bar(y)),
        // using the definition of the mean, and simplifying
        cov_xy = (cov_xy/static_cast<double> (length)) - mean_x*mean_y;
        b1_ = cov_xy/var_x;
        b0_  = mean_y - b1_*mean_x;
        return 0; 
    }
    /// Computes the linear regression coefficients
    int computeLinearRegressionCoeffs(const int length, const float pSrcY[])
    {
        double cov_xy;
        double mean_x;
        double mean_y;
        double var_x;
        float mean_y32;
        // Mean of x - analytic formula for evenly spaced samples starting
        // at indx 0.  This is computed by simplifying Gauss's formula.
        auto len64 = static_cast<uint64_t> (length);
        mean_x = 0.5*static_cast<double> (len64 - 1); 
        // Note, the numerator is the sum of consecutive squared numbers.
        // In addition we simplify.
        var_x = static_cast<double> ( ((len64 - 1))*(2*(len64 - 1) + 1) )/6.
              - mean_x*mean_x;
        ippsMean_32f(pSrcY, length, &mean_y32, ippAlgHintAccurate); 
        mean_y = static_cast<double> (mean_y32);
        cov_xy = 0;
        #pragma omp simd reduction(+:cov_xy)
        for (int i=0; i<length; i++)
        {
            cov_xy = cov_xy
                   + static_cast<double> (i)*static_cast<double> (pSrcY[i]);
        }
        // This is computed by expanding (x_i - bar(x))*(y_i - bar(y)),
        // using the definition of the mean, and simplifying
        cov_xy = (cov_xy/static_cast<double> (length)) - mean_x*mean_y;
        b1_ = cov_xy/var_x;
        b0_  = mean_y - b1_*mean_x;
        return 0; 
    }
    /// Sets the intercept 
    void setIntercept(const double b0)
    {
        b0_ = b0;
    }
    /// Sets the slope
    void setSlope(const double b1)
    {
        b1_ = b1;
    }
//private:
    /// Detrend parameters
    DetrendParameters parms_;
    /// The y-intercept
    double b0_ = 0;
    /// The slope
    double b1_ = 0;
    /// Flag indicating this is initialized
    bool linit_ = true;
};
//============================================================================//

DetrendParameters::DetrendParameters(const RTSeis::Precision precision) :
    pImpl(std::make_unique<DetrendParmsImpl>())
{
    pImpl->precision_ = precision;
}

DetrendParameters::DetrendParameters(const DetrendParameters &parameters)
{
    *this = parameters;
}

DetrendParameters&
    DetrendParameters::operator=(const DetrendParameters &parameters)
{
    if (&parameters == this){return *this;}
    //if (pImpl){pImpl->clear();}
    //pImpl = std::make_unique<DetrendParms>(*parameters.pImpl);
    pImpl = std::unique_ptr<DetrendParmsImpl>
                         (new DetrendParmsImpl(*parameters.pImpl));
    return *this;
}

DetrendParameters::~DetrendParameters() = default;

void DetrendParameters::clear()
{
    pImpl->clear();
}

RTSeis::Precision DetrendParameters::getPrecision() const
{
    return pImpl->precision_;
}

RTSeis::ProcessingMode DetrendParameters::getProcessingMode() const
{
    return pImpl->mode_;
}

bool DetrendParameters::isInitialized() const
{
    return pImpl->linit_;
}
//============================================================================//
Detrend::Detrend() :
    pDetrend_(std::make_unique<DetrendImpl>())
{
    clear();
}

Detrend::Detrend(const Detrend &detrend)
{
    *this = detrend;
}

Detrend::Detrend(const DetrendParameters &parameters) : 
    pDetrend_(std::make_unique<DetrendImpl>())
{
    setParameters(parameters);
}

Detrend& Detrend::operator=(const Detrend &detrend)
{
    if (&detrend == this){return *this;}
    if (pDetrend_){pDetrend_->clear();}
    //pDetrend_ = std::make_unique<DetrendImpl> (*detrend.pDetrend_);
    pDetrend_ = std::make_unique<DetrendImpl> (*detrend.pDetrend_);
    return *this;
}

Detrend::~Detrend()
{
    clear();
}

void Detrend::clear()
{
    pDetrend_->clear();
}

void Detrend::setParameters(const DetrendParameters &parameters)
{
    pDetrend_->clear();
    if (!parameters.isInitialized())
    {
        RTSEIS_THROW_IA("%s", "Detrend parameters are malformed");
    }
    pDetrend_->setParameters(parameters);
}

void Detrend::apply(const int nx, const double x[], double y[])
{
    pDetrend_->setIntercept(0);
    pDetrend_->setSlope(0);
    if (nx <= 0){return;} // Nothing to do
    if (nx < 2 || x == nullptr || y == nullptr)
    {
        if (nx < 2){RTSEIS_THROW_IA("%s", "At least 2 points required");}
        if (x == nullptr){RTSEIS_THROW_IA("%s", "x is null");}
        RTSEIS_THROW_IA("%s", "y is null");
    }
    pDetrend_->apply(nx, x, y);
}

void Detrend::apply(const int nx, const float x[], float y[])
{
    pDetrend_->setIntercept(0);
    pDetrend_->setSlope(0);
    if (nx <= 0){return;} // Nothing to do
    if (nx < 2 || x == nullptr || y == nullptr)
    {
        if (nx < 2){RTSEIS_THROW_IA("%s", "At least 2 points required");}
        if (x == nullptr){RTSEIS_THROW_IA("%s", "x is null");}
        RTSEIS_THROW_IA("%s", "y is null");
    }
    pDetrend_->apply(nx, x, y);
}
